///usage:  ini ROOT528
///usage:  root -l -q  readtribn-data.root 'fitTemplates.C(0,0,3)' ## Only bkg models+data
///usage:  root -l -q  readtribn-data.root readtribn-signal.root 'fitTemplates.C(0,0,3)' ## bkg+signal+data
///usage:  root -l -q  readtribn-data.root readtribn-signal.root 'fitTemplates.C(0,0,3,"merge_templates.txt")' # if you want merge sometemplates 
///usage:	remember !!!! -> first root ntuple must be a ntuple with data and bkg templates !!!!!
///usage:
///usage:	fitTemplates.C(_novl,_corr,_scenario)
///usage:         Parameters of the function:
///usage:         _novl = 0,1 (TCHPT, TCHP6), Btag operation point of mass templates
///usage:         _corr= 0,1 (uncorrected,corrected) bb-purity correction application
///usage:         _scenario = 0,1,2,3
///usage:                 * 0 -- Standard unbinned ML
///usage:                 * 1 -- Standard  binned Chi2
///usage:                 * 2 -- Unbinned ML with help of RooMinuit
///usage:                 * 3 -- Binned Chi2 with help of RooMinuit

///usage:
///usage:

///usage:  We have 9 bkg templates:
///usage:	bbB bBb Bbb
///usage:	bbC bCb Cbb
///usage:	bbQ bQb Qbb
///usage:	+ 1 signal template
///usage:  strategy of combining:
///usage:  step1: Bbb+bBb
///usage:  step1: Cbb+bCb
///usage:  step1: Qbb+bQb

///usage:
///usage:
///usage:   if correlated
///usage:  step2: bbQ+bbC


///usage: script can write all output to log file 
///usage: you can extract intresing points by doing
///usage: sed -ne '/\+MYOUTPUT/,/\-MYOUTPUT/p' ROOFIT.log 




#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooConstVar.h"
#include "RooProdPdf.h"
#include "RooAddPdf.h"
#include "RooMinuit.h"
#include "RooFitResult.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "TH1.h"
#include "RooHistPdf.h"
#include "RooPlot.h"
#include "RooNLLVar.h"
#include "RooChi2Var.h"

#include "TObjString.h"
#include "TList.h"
#include "TMath.h"
#include "TLegend.h"
#include "TGraph.h"
#include "TArrayD.h"

using namespace RooFit ;

#include <iomanip>
#include <string>

#include <unistd.h>
#include <stdio.h>


#include "TGraphAnalysis.h"

//#define doLogFile ///usage
#define LOGFILE "ROOFIT.log" ///usage


///usage: to work with switches:
///usage: sed -e 's/#define doMIGRAD/\/\/\#define doMIGRAD/g' fitTemplates.C > myFit.C.tmp; mv myFit.C.tmp fitTemplates.C
///usage:  sed -e 's/\/\/\#define doHESSE/#define doHESSE/g' fitTemplates.C > myFit.C.tmp; mv myFit.C.tmp fitTemplates.C;

#define doMIGRAD ///usage
#define doHESSE ///usage 
#define doMINOS ///usage 

///usage: if you want to combine all relevevant info in the tables;
///usage: sed -ne '/\+MYOUTPUT/,/\-MYOUTPUT/p' ROOFIT.log | grep -E "=|VALS|MIGRAD|MINOS|RATIOS" | tr  '=' ';' | tr -d ' ' | sed -e 's/_TCHPTTrig0//g' >| table.csv 
///usage: Make_table "my table" table.csv ascii
///usage:
///usage:
///usage:


///To merge histograms
#include "templateMerger.h"

int colors [] =
{
	kRed,
	kGreen,
	kMagenta,
	kYellow,
	kRed-2,
	kGreen-2,
	kMagenta-2,
	kYellow-2,
	kRed+3,
	kGreen+3,
	kMagenta+3,
	kYellow+3,
	

}

///Templates && data histogram properties
TH1 * dataHist=0;
double minVal=-1e30;
double maxVal=+1e30;
double dataYield=1e0;
int bins = 0;

///Flag to understand the status of the current code execution
bool _isOk=false;



////Containers

TList * arrayOfHistos=0; /// First element is always data histogram
TList * arrayOfNames =0; 
TList * arrayOfRooDataHist=0;
TList * arrayOfRooHistPdf=0;
TList * arrayOfCanvas=0; 
TList * arrayOfYields=0;
TList * arrayOfLegends=0;


///Instead of direct histogram names

//const  char * histonames[] = {
//  "mjjtchp6Trig0", //data
//  "mjjtchp6TempCat0UdsTrig0",  // Qbb
//  "mjjtchp6TempCat2UdsTrig0"//,  // bbQ
//  "mjjtchp6TempCat0CTrig0",    // Cbb
//  "mjjtchp6TempCat0BTrig0",    // Bbb
//  "mjjtchp6TempCat2BTrig0"    // bbB
//};



/// we will use two additional index arrays (with the same size) carring info about
/// 1) flavor 0,1,2 (for "Q", "C", "B");  size of values == nfc 
/// 2) Cat 0,1,2 ( for Bbb,bBb,bbB); size of values == ncateg
/// sizes of flavors and positions arrays must be the same

/// Here   


///	Bbb bBb bbb
///	Cbb bCb bbC
///	Qbb bQb bbQ

/// is coded

/// Carries 'Flavor' index
 const int flavors [] = {
	2, //B
//	2, //B
//	2, //B
//	1, //C
	1, //C
//	1,  //C
//	0, //Q
	0, //Q
//	0  //Q
	};

/// Carries 'Cat' index
 const int positions [] = {
	0, //Bbb
//	1, //bBb
//	2, //bbB
//	0, //Cbb
	1, //bCb
//	2, //bbC
//	0, //Qbb
	1, //bQb
//	2  //bbQ
	};

/**

		Chi2 Test function
		my own code
	
		Parameters of the function

		_roohistpdfs = TList* of TH1* of templates. First ellement must be data histo
		_yelds =  TArrayD* of yields
		

	
**/



double   chi2test2(TList * _roohistpdfs, TArrayD * _yelds)
{

double res=-1e10;
if (!(_roohistpdfs && _yelds)!=0) {std::cout<<"chi2test :: empty pointers, error"<<std::endl; return res;}
if ( _roohistpdfs->GetSize()<2 ) { std::cout<<"chi2test :: expect to have data th1 and at least 1 MC template, error"<<std::endl; return res;}


TH1 * _dataHist = (TH1*) _roohistpdfs->At(0);
TH1F * _combinedHist = new TH1F("combined_hist","combined_hist",_dataHist->GetNbinsX(),_dataHist->GetXaxis()->GetXmin(),_dataHist->GetXaxis()->GetXmax());

cout<<"chi2Test validation: "<<endl;
cout<<"Data integral "<<_dataHist->Integral()<<endl;

//_dataHist->Sumw2();
//_combinedHist->Sumw2();


double sum=0;
double _err[500];
int degree=_dataHist->GetNbinsX();

for (int kk=1;kk<=_dataHist->GetNbinsX();kk++) _err[kk] =0;

for (int k=1;k<_roohistpdfs->GetSize();k++) {
degree--;
/*cout<<"I'm " << ((TH1 * ) _roohistpdfs->At(k))->GetName()<<" with integral "<<((TH1 * ) _roohistpdfs->At(k))->Integral() <<" and yield "<<_yelds->GetAt(k)<<endl;*/

_combinedHist->Add((TH1 * ) _roohistpdfs->At(k),_yelds->GetAt(k));
//for (int kk=1;kk<=_dataHist->GetNbinsX();kk++) {
/*
cout<<"kk = "<<kk<<endl;
cout<<"case1 = "<<pow(_yelds->GetAt(k),2)*pow(((TH1 * )_roohistpdfs->At(k))->GetBinError(kk),2)<<endl;
cout<<"case2 = "<<pow(_yelds->GetAt(k)*((TH1 * )_roohistpdfs->At(k))->GetBinError(kk),2)<<endl;
cout<<"case3 = "<<_yelds->GetAt(k)*((TH1 * )_roohistpdfs->At(k))->GetBinError(kk)<<endl;
*/
//_err[kk]+=pow(_yelds->GetAt(k),2)*pow(((TH1 * )_roohistpdfs->At(k))->GetBinError(kk),2);
//cout<< "Bin Erorr = "<< ((TH1 * )_roohistpdfs->At(k))->GetBinError(kk)<<" yield = "<< _yelds->GetAt(k) <<" _err = "<<_err[kk]<<endl;
//}
sum+=_yelds->GetAt(k);
}


for (int kk=1;kk<=_dataHist->GetNbinsX();kk++) {
/*
cout<<"kk = "<<kk<<endl;
cout<<"case1 = "<<pow(_yelds->GetAt(k),2)*pow(((TH1 * )_roohistpdfs->At(k))->GetBinError(kk),2)<<endl;
cout<<"case2 = "<<pow(_yelds->GetAt(k)*((TH1 * )_roohistpdfs->At(k))->GetBinError(kk),2)<<endl;
cout<<"case3 = "<<_yelds->GetAt(k)*((TH1 * )_roohistpdfs->At(k))->GetBinError(kk)<<endl;
*/
//_err[kk]+=pow(_combinedHist->GetBinError(kk),2);
_err[kk]+=pow(_combinedHist->GetBinError(kk),2);
//cout<<"Bin = "<<kk<<endl;
//cout<< "Bin Erorr^2 = "<<_err[kk]<<endl;
//cout<< "Bin Erorr = "<<_combinedHist->GetBinError(kk)<<endl;
//cout<< "Bin Content = "<<_combinedHist->GetBinContent(kk)<<endl;
//cout<< "Bin Erorr = "<< ((TH1 * )_roohistpdfs->At(k))->GetBinError(kk)<<" yield = "<< _yelds->GetAt(k) <<" _err = "<<_err[kk]<<endl;
}
 

cout<<"Total MC events is "<<sum<<endl;
if (_dataHist !=0 ) {

//cout<<"I'm combined histo  with integral "<<_combinedHist->Integral()<<endl;

//_dataHist->Sumw2();
//_combinedHist->Sumw2();

double _chi2=0e0;
double _diff=0e0;

for (int kk=1;kk<=_dataHist->GetNbinsX();kk++)
{

//double _error =  _err[kk];
double _error =  _dataHist->GetBinContent(kk) + _err[kk];
//double _error =  _dataHist->GetBinContent(kk) + _combinedHist->GetBinContent(kk);


_diff = _dataHist->GetBinContent(kk) - _combinedHist->GetBinContent(kk);
double _add = _diff*_diff/_error;

/*
std::cout<<"I'm bin "<<kk<<" with data "<<_dataHist->GetBinContent(kk)<<" and MC "<<_combinedHist->GetBinContent(kk)<<std::endl;
std::cout<<" and diff "<<_diff <<" and add "<<_add<<std::endl;	
*/
	 if (_combinedHist->GetBinContent(kk)>0) _chi2+=_add;
	else	
	if (_dataHist->GetBinContent(kk)>0)  _chi2+=_diff*_diff/_dataHist->GetBinContent(kk);


if (_combinedHist->GetBinContent(kk)==0 ) degree--;
}

}

cout<<"4.1) My Chi2 = "<<_chi2<<endl;
cout<<"4.2) My ndof = "<<degree<<endl;
cout<<"4.2) My Chi2/ndof = "<<_chi2/degree<<endl;

delete _combinedHist;
return _chi2;


}

/**	
		Doesn't work correctly (:-
	
		It's kept to remember some ideas of the use TH1::Chi2Test()
		
		Chi2 Test function
		based on TH1::Chi2Test()
	
		Parameters of the function

		_roohistpdfs = TList* of TH1* of templates. First ellement must be data histo
		_yelds =  TArrayD* of yields
		_options = TString of options	

	
**/


double   chi2test(TList * _roohistpdfs, TArrayD * _yelds, TString _options)
{

double res=-1e10;
if (!(_roohistpdfs && _yelds)!=0) {std::cout<<"chi2test :: empty pointers, error"<<std::endl; return res;}
if ( _roohistpdfs->GetSize()<2 ) { std::cout<<"chi2test :: expect to have data th1 and at least 1 MC template, error"<<std::endl; return res;}


TH1 * _dataHist = (TH1*) _roohistpdfs->At(0);
TH1F * _combinedHist = new TH1F("combined_hist","combined_hist",_dataHist->GetNbinsX(),_dataHist->GetXaxis()->GetXmin(),_dataHist->GetXaxis()->GetXmax());

cout<<"chi2Test validation: "<<endl;
cout<<"Data integral "<<_dataHist->Integral()<<endl;

 
double sum=0;
for (int k=1;k<_roohistpdfs->GetSize();k++) {
cout<<"I'm " << ((TH1 * ) _roohistpdfs->At(k))->GetName()<<" with integral "<<((TH1 * ) _roohistpdfs->At(k))->Integral() <<" and yield "<<_yelds->GetAt(k)<<endl;
_combinedHist->Add((TH1 * ) _roohistpdfs->At(k),_yelds->GetAt(k));
sum+=_yelds->GetAt(k);
}

cout<<"Total MC events is "<<sum<<endl;

for (int kk=1;kk<=_dataHist->GetNbinsX();kk++)
if (_combinedHist->GetBinContent(kk)==0 && _dataHist->GetBinContent(kk)==0) {
	_combinedHist->SetBinContent(kk,1);
	_dataHist->SetBinContent(kk,1);
}


if (_dataHist !=0 ) {

cout<<"I'm combined histo  with integral "<<_combinedHist->Integral()<<endl;

_dataHist->Sumw2();
_combinedHist->Sumw2();

 if (_options.Length() == 0 ) res=_dataHist->Chi2Test( _combinedHist,"NORM OF UF P CHI2");
else  res=_dataHist->Chi2Test( _combinedHist,_options.Data());

}

delete _combinedHist;
return res;
}



/**
	Parameters of the function:
	_novl = 0,1 (TCHPT, TCHP6), Btag operation point of mass templates
	_corr= 0,1 (uncorrected,corrected) bb-purity correction application


	Please verify that: gROOT->GetListOfFiles()->GetSize()>1 otherwise signal template is not added

	

**/


bool selectSignalTemplates(int _novl=0)
{

std::string tSel("Trig1");

if (!arrayOfHistos) 
{arrayOfHistos = new TList();
arrayOfHistos->SetOwner(kTRUE); ///To auto-delete all content of the container
}
if (!arrayOfNames) {
arrayOfNames = new TList();
arrayOfNames->SetOwner(kTRUE); ///To auto-delete all content of the container
}

///Read histograms from file	
TFile * file1=0;
 if (gROOT->GetListOfFiles()->GetSize()>1) file1=(TFile*) gROOT->GetListOfFiles()->At(1);

///Parameters of reading 
  const int novl = 2;

  TH1F* mjjTemplate[novl];
  const std::string sovl[novl] = { "TCHPT", "TCHP6" };

///Reading histo MC (templates) for signal

_isOk=true; /// assume that all histograms are available
   
for (int iovl=0; iovl<novl; ++iovl)    
{
///mjjbtTCHPTTrig1
if (file1!=0) mjjTemplate[iovl] = (TH1F*) file1->Get(Form("mjjbt%s%s",sovl[iovl].c_str(),tSel.c_str()));
else
{
 mjjTemplate[iovl]  =0 ;

}


}


if (mjjTemplate[_novl] !=0)
	{
 		mjjTemplate[_novl]->SetDirectory(gROOT); ///To make histograms independent on _file1;
		arrayOfHistos->Add(mjjTemplate[_novl]);


///Form Name of the template:		
	 arrayOfNames->Add(new TObjString(Form("signal_%s%s",sovl[_novl].c_str(),tSel.c_str())));
	}
else {
std::cout << "Histogram not found: " << Form("mjjbt%s%s",sovl[_novl].c_str(),tSel.c_str())<< std::endl;
_isOk=false;
}


if (!_isOk ) {std::cout<<"There are problems to read all needed signal templates from the file, error"<<std::endl;
return false;}
_isOk=false; /// reset to ini value


return true;
}



/**
	Parameters of the function:
	_novl = 0,1 (TCHPT, TCHP6), Btag operation point of mass templates
	_corr= 0,1 (uncorrected,corrected) bb-purity correction application


**/


bool selectTemplates(int _novl=0, int _corr=0)
{

std::string tSel("Trig0");
if (!arrayOfHistos) {
arrayOfHistos = new TList();
arrayOfHistos->SetOwner(kTRUE); ///To auto-delete all content of the container
}
if (!arrayOfNames) {
arrayOfNames = new TList();
arrayOfNames->SetOwner(kTRUE); ///To auto-delete all content of the container
}
///Read histograms from file	
TFile * file0=0;
 if (gROOT->GetListOfFiles()->GetSize()>0) file0=(TFile*) gROOT->GetListOfFiles()->At(0);


///Parameters of reading 
  const int ncateg = 3;
  const int ncorr=2;
  const int ntpat=1;
  const int nfc = 3;
  const int novl = 2;

  TH1F* mjjTemplate[nfc][novl][ncateg][ncorr][ntpat];
//  const std::string tFlav[nfc] = {"Uds", "C", "B"};
  const std::string sFlav[nfc] = {"Q", "C", "B"};

  const std::string sfc[nfc] = { "q", "c", "b" };
  const std::string sovl[novl] = { "TCHPT", "TCHP6" };

///Reading Data

if (file0==0) { cout<<"Can't read data, error"<<endl;return false;}

 TH1F* mjjdata[novl];
  for (int iovl=0; iovl<novl; ++iovl) 
  if (file0!=0)  mjjdata[iovl] = (TH1F*) file0->Get(Form("mjjbt%s%s",sovl[iovl].c_str(),tSel.c_str()));
  else mjjdata[iovl] =0;

///Putting Data into arrayOfHistos
for (int iovl=0; iovl<novl; ++iovl) 
if (mjjdata[iovl]!=0)
{
	mjjdata[iovl]->SetDirectory(gROOT); ///To make histograms independent on _file0;
	if (iovl==_novl)  {
		arrayOfHistos->Add(mjjdata[iovl]);	
		arrayOfNames->Add(new TObjString(Form("data_%s%s",sovl[iovl].c_str(),tSel.c_str())));
		_isOk=true;
	}

}

if (!_isOk) {std::cout<<"The data mass spectrum of "<<sovl[_novl]<<" operation point was not found, error"<<std::endl;
return false;
}
_isOk=false; ///reset to ini value


///Reading histo MC (templates)

_isOk=true; /// assume that all histograms are available
  for (int ifc=0; ifc<nfc; ++ifc) 
    for (int iovl=0; iovl<novl; ++iovl) 
      for (int icateg=0; icateg<ncateg; ++icateg) 
        for (int icorr=0; icorr<ncorr; ++icorr) 
          for (int itpat=0; itpat<ntpat; ++itpat) 
{

if (file0!=0) mjjTemplate[ifc][iovl][icateg][icorr][itpat] = (TH1F*) file0->Get(Form("mjjTemp_%s_%s_Cat%dCorr%dTpat%d%s",sfc[ifc].c_str(),sovl[iovl].c_str(),icateg,icorr,itpat,tSel.c_str()));
else
{
 mjjTemplate[ifc][iovl][icateg][icorr][itpat]  =0 ;

}


}



for (unsigned int i=0;i<sizeof(flavors)/sizeof(int);i++) {  
	if (mjjTemplate[flavors[i]][_novl][positions[i]][_corr][0]!=0)
	{
	 mjjTemplate[flavors[i]][_novl][positions[i]][_corr][0]->SetDirectory(gROOT); ///To make histograms independent on _file0;		
	arrayOfHistos->Add(mjjTemplate[flavors[i]][_novl][positions[i]][_corr][0]);

///Form Name of the template:
	TString _tmpl = "bbb";
	_tmpl.Replace(positions[i],1,sFlav[flavors[i]].c_str());
	
	 arrayOfNames->Add(new TObjString(Form("%s_%s%s",_tmpl.Data(),sovl[_novl].c_str(),tSel.c_str())));	

	} 

	else {
 std::cout << "Histogram not found: " << Form("mjjTemp_%s_%s_Cat%dCorr%dTpat%d%s",sfc[flavors[i]].c_str(),sovl[_novl].c_str(),positions[i],_corr,0,tSel.c_str())<< std::endl;
_isOk=false;
	}

}
if (!_isOk ) {std::cout<<"There are problems to read all needed templates from the file, error"<<std::endl;
return false;}
_isOk=false; /// reset to ini value


return true;
}

/**
		Performs merging 
		_mother == mother template ( ie bbB)
		_daughters == TList * of daughter templates (TObjString* names)
		_mode == mode of operation: Merge && Normalize -- 1, Normalize && Merge && Normalize 
	It changes :
		* arrayOfHistos
		* arrayOfNames
		deletes  corresponded _daughters templates references
	

**/

bool templateMerge(TString _mother, TList *_daughters, int mode=1)
{

	if ( _mother.Length()==0 && _daughters == 0 && _daughters->IsEmpty()) {std::cout<<"templateMerge --> something wrong"<<std::endl;return false;}

TList *_remove_arrayOfHistos = new TList();
TList *_remove_arrayOfNames = new TList();


///Find mother
TH1 *  _motherTh1=0;
for (int i=0; i<arrayOfNames->GetSize();i++)
{
TString _name = ( (TObjString *) arrayOfNames->At(i))->GetString();
if (_name.Contains(_mother)){ _motherTh1 =  (TH1*)arrayOfHistos->At(i); break;}
}


if (_motherTh1==0) {std::cout<<"templateMerge --> something wrong mother"<<std::endl;return false;}

///Find daughter templates
for (int i=0; i<arrayOfNames->GetSize();i++)
{

TH1 * tempHis = (TH1*)arrayOfHistos->At(i);

TString _name = ( (TObjString *) arrayOfNames->At(i))->GetString();
 for (int jj=0; jj<_daughters->GetSize();jj++)
{
	TString _nameDaughter =  ( (TObjString *) _daughters->At(jj))->GetString();
	if (_name.Contains(_nameDaughter)) {

	_remove_arrayOfNames->Add(arrayOfNames->At(i));
	_remove_arrayOfHistos->Add(arrayOfHistos->At(i));

	} ///if
} ///for


} ///for


///Merging daughters to mother && Deleting daughter from arrayOfNames arrayOfHistos
if (!templateMerger(_motherTh1,_remove_arrayOfHistos,mode)) {std::cout<<"templateMerge --> something wrong with merging"<<std::endl;return false;}


for (int i=0; i<_remove_arrayOfNames->GetSize();i++)
{
	arrayOfHistos->Remove( _remove_arrayOfHistos->At(i));
	delete  _remove_arrayOfHistos->At(i);
	arrayOfNames->Remove( _remove_arrayOfNames->At(i));
	delete _remove_arrayOfNames->At(i);
}


delete _remove_arrayOfHistos;
delete _remove_arrayOfNames;
return true;

}




/**
        Parameters of the function:
        _novl = 0,1 (TCHPT, TCHP6), Btag operation point of mass templates
        _corr= 0,1 (uncorrected,corrected) bb-purity correction application
	_scenario = 0,1,2,3
		* 0 -- Standard unbinned ML
		* 1 -- Standard  binned Chi2
		* 2 -- Unbinned ML with help of RooMinuit
		* 3 -- Binned Chi2 with help of RooMinuit

	_merging=""
		* name of the file having info about merging (possible content):

		#comments here
		#mother daughter1 daughter2 daughter3 etc
		Bbb	bBb	bbB	
		Qbb	bQb	bbQ
		#etc
**/


void fitTemplates(int _novl=0, int _corr=0, int _scenario =3, TString _merging="") ///usage
{
#ifdef doLogFile
freopen(LOGFILE,"w",stdout);
dup2(fileno(stdout), fileno(stderr));
#endif

///Read all needed stuff from the _file0 
if (!selectTemplates(_novl,_corr) ) {std::cout<<"Something wrong with reading bkg and data, return"<<std::endl; return;}

if (!selectSignalTemplates(_novl) ) {std::cout<<"Something wrong with reading signal"<<std::endl; }

if (arrayOfHistos && !arrayOfHistos->IsEmpty()) dataHist=(TH1*)arrayOfHistos->At(0);



if (dataHist!=0)
{
///To solve the problem zero statistics in bin
for (int kk=1;kk<=bins;kk++)  if (dataHist->GetBinError(kk)==0 && dataHist->GetBinContent(kk)==0) dataHist->SetBinError(kk,1);
	minVal=dataHist->GetXaxis()->GetXmin();
	maxVal=dataHist->GetXaxis()->GetXmax();
	bins = dataHist->GetNbinsX();
	dataYield=dataHist->Integral();

	///Validation
	cout<<"Validation of mass data spectrum: "<<endl;
	cout<<"MinVal = "<<minVal<<endl;
	cout<<"MaxVal = "<<maxVal<<endl;
	cout<<"Number of bins = "<<bins<<endl;
	 cout<<"Number of all events (data) = "<<dataYield<<endl;

/*double sum=0;
for (int kk=1;kk<=bins;kk++) {
if (dataHist->GetBinError(kk)==0 && dataHist->GetBinContent(kk)==0)dataHist->SetBinError(kk,1);
cout<<"Bin # "<<kk<<" val ="<<dataHist->GetBinContent(kk)<<" error= "<<dataHist->GetBinError(kk)<<endl;
sum+=dataHist->GetBinContent(kk);

}

cout<<"sum = "<<sum<<endl;
*/
} else {cout<<"dataHist is empty, return"<<endl; return;}



RooRealVar Mjj("Mjj","Mjj",minVal,maxVal) ;
Mjj.setBins(bins);


arrayOfRooDataHist=new TList();
arrayOfRooDataHist->SetOwner(kTRUE); ///To auto-delete all content of the container

arrayOfRooHistPdf=new TList();
arrayOfRooHistPdf->SetOwner(kTRUE); ///To auto-delete all content of the container

arrayOfCanvas=new TList();
arrayOfCanvas->SetOwner(kTRUE); ///To auto-delete all content of the container

arrayOfYields = new TList();
arrayOfYields->SetOwner(kTRUE); ///To auto-delete all content of the container

arrayOfLegends=new TList();
arrayOfLegends->SetOwner(kTRUE); ///To auto-delete all content of the container



///Here possible merging

if (_merging.Length()>0) {

ifstream f(_merging);

///Parsing of file:
TString  str;

///Read first line
str.ReadLine(f);

while (str.Length()>0 )
{
TObjArray* tokens=str.Tokenize(TString(' '));

if (tokens->GetEntries()>1 && !str.Contains("#")) {



TString mother = ((TObjString*) tokens->At(0))->GetString();
TList * daughters = new TList();
daughters->SetOwner(kTRUE); ///To auto-delete all content of the container

for (int dd=1;dd<tokens->GetEntries();dd++) 	daughters->Add(new TObjString(tokens->At(dd)->GetName()));
	

if ( !templateMerge(mother, daughters, 0)) cout<<"Something wrong with merging!!!"<<endl;
//if ( !templateMerge(mother, daughters, 1)) cout<<"Something wrong with merging!!!"<<endl;

///old code
///Case Bbb+bbB+Cbb+bbQ->Bbb
/*
mother = "Bbb";
daughters->Add(new TObjString("bbB"));
daughters->Add(new TObjString("Cbb"));
daughters->Add(new TObjString("bbQ"));
if ( !templateMerge(mother, daughters, 0)) cout<<"Something wrong with merging!!!"<<endl;
daughters->Clear();
*/

///End of merging
daughters->Clear();
delete daughters;
delete tokens;
} /// token
str.ReadLine(f); 
} 
} ///if _merging


//for (unsigned int i=0;i<sizeof(histonames)/sizeof(char *);i++)
for (int i=0; i<arrayOfNames->GetSize();i++)
{
TH1 * tempHis = (TH1*)arrayOfHistos->At(i);
TString _name = ( (TObjString *) arrayOfNames->At(i))->GetString();



if (tempHis!=0) 
{	
/*
TCanvas * can =new TCanvas(); 
tempHis->Draw();
*/
arrayOfRooDataHist->Add(new  RooDataHist(_name+TString("_histdata"),_name+TString("_histdata"),Mjj,tempHis));
arrayOfRooHistPdf->Add( new RooHistPdf(_name+TString("_histpdf"),_name+TString("_histpdf"),RooArgList(Mjj),*((RooDataHist*)arrayOfRooDataHist->Last()),0)); ///Last argument is a degree of internal interpolation
//cout<<_name<<endl;
//if (_name.Contains("Bbb_")) arrayOfYields->Add( new RooRealVar(_name+TString("_yield"),_name+TString("_yield"),3331e0,0.,dataYield) );
//else if (_name.Contains("bbB_")) arrayOfYields->Add( new RooRealVar(_name+TString("_yield"),_name+TString("_yield"),7608e0,0.,dataYield) );
//else arrayOfYields->Add( new RooRealVar(_name+TString("_yield"),_name+TString("_yield"),dataYield/5e0,0.,dataYield) );
arrayOfYields->Add( new RooRealVar(_name+TString("_yield"),_name+TString("_yield"),dataYield/3e0,0.,dataYield) );
//arrayOfYields->Add( new RooRealVar(_name+TString("_yield"),_name+TString("_yield"),dataYield/5e0,0.,dataYield) );
//arrayOfYields->Add( new RooRealVar(_name+TString("_yield"),_name+TString("_yield"),dataYield,0e0,1e20,"") );
//arrayOfYields->Add( new RooRealVar(_name+TString("_yield"),_name+TString("_yield"),0e0,0e0,dataYield*20,"") );
//arrayOfYields->Add( new RooRealVar(_name+TString("_yield"),_name+TString("_yield"),1e2,0e0,1e30,"") );
arrayOfCanvas->Add (new  TCanvas(_name+TString("_canva"),_name+TString("_canva")) ); /// skip canva for data


arrayOfLegends->Add(new TLegend(0.4,0.6,0.89,0.89));


}

}



///For plotting combined templates
arrayOfCanvas->Add (new TCanvas(TString("combined")+TString("_canva"),TString("combined")+TString("_canva")));
arrayOfLegends->Add(new TLegend(0.4,0.6,0.89,0.89));
TGraph _combined;
TList * _combinedList = new TList();
_combinedList->SetOwner(kTRUE); 


///For plotting final result
arrayOfCanvas->Add (new TCanvas(TString("final")+TString("_canva"),TString("final")+TString("_canva")));
arrayOfLegends->Add(new TLegend(0.4,0.6,0.89,0.89));

 
///Validation of histPdf
cout<<"Validation: "<<endl;
//for (unsigned int i=0;i<sizeof(histonames)/sizeof(char *);i++)
for (int i=0; i<arrayOfNames->GetSize();i++)
{

if (i==0) continue; /// skip data

TString _name = ( (TObjString *) arrayOfNames->At(i))->GetString();
((TCanvas*)arrayOfCanvas->At(i))->cd();
TLegend *  leg = (TLegend*) arrayOfLegends->At(i);


TH1F* histPdf=(TH1F*) ((RooHistPdf*)arrayOfRooHistPdf->At(i))->createHistogram(_name+TString("_histpdf_validation"),Mjj));
if (histPdf!=0) {
	cout<<"Integral is "<<histPdf->Integral("")<<endl; /// or ->Integral("width")?
	histPdf->SetLineColor(colors[i]);
	leg->AddEntry(histPdf,_name,"l");
	histPdf->Draw();
	leg->Draw();
	_combined.SetPoint(i,i,histPdf->GetMaximum());	
//	cout<<"name  "<<histPdf->GetName()<<endl;
	_combinedList->Add(histPdf);
	}


}

///Validation of combined hist
/*
for (int i=0; i< _combinedList->GetSize();i++)
{
TH1F* histPdf=(TH1F*) _combinedList->At(i);
 cout<<"I'm histo "<<histPdf->GetName()<<endl;

double sum=0;
if (histPdf) for (int kk=1;kk<=bins;kk++) {
cout<<"Bin # "<<kk<<" val ="<<histPdf->GetBinContent(kk)<<" error= "<<histPdf->GetBinError(kk)<<endl;
sum+=histPdf->GetBinContent(kk);
}

cout<<"sum = "<<sum<<endl;

}
*/
///Plot combined 
((TCanvas*)arrayOfCanvas->FindObject(TString("combined")+TString("_canva")))->cd();
///Make Sorting
_combined.Sort(&TGraph::CompareY,kFALSE); ///descending on Y

Double_t * templPos = _combined.GetX();
_isOk=true; ///flag the histogram with the highest GetMaximum()

for (int i=0; i<_combined.GetN();i++)
{
if (templPos[i] == 0 ) continue; /// skip possible data

TLegend *  leg = (TLegend*) arrayOfLegends->At(arrayOfLegends->GetSize()-2); ///prelast;
TString _name = ( (TObjString *) arrayOfNames->At(templPos[i]))->GetString();


	TH1F* histPdf=(TH1F*) _combinedList->FindObject(_name+TString("_histpdf_validation__Mjj"));

	if (histPdf!=0) {
		 leg->AddEntry(histPdf,_name,"l");
		histPdf->SetLineColor(colors[(Int_t)templPos[i]]);
	if (_isOk) {
		histPdf->Draw();
		_isOk=false;	

		}
	else 		histPdf->Draw("same");
		
		leg->Draw();

	
	}


}


///Create Model Pdf
// ---------------------------------------------------------------------

RooArgList * yields = new RooArgList();
RooArgList * shapes = new RooArgList();

//for (unsigned int i=0;i<sizeof(histonames)/sizeof(char *);i++)
for (int i=0; i<arrayOfNames->GetSize();i++)
{

if (i==0)  continue; /// skip data

cout<<"Adding "<<((RooRealVar*)arrayOfYields->At(i))->GetName()<<endl;

yields->add(*(RooRealVar*)arrayOfYields->At(i));
shapes->add(*(RooHistPdf*)arrayOfRooHistPdf->At(i));

}

// Sum all components into a composite model p.d.f.
// ---------------------------------------------------------------------

 RooAddPdf  model("model","model",*shapes,*yields);

RooFitResult*  result=0;



//  F i t   a n d   p l o t   e x t e n d e d   m o d e l and LatexTables
// ---------------------------------------------------------------------

// Make list of model parameters
  RooArgSet* params = model.getParameters(Mjj) ;

//METHOD 1 : Standard unbinned ML
// ---------------------------------------------------------------------

if (_scenario == 0)
{


//RooFitResult*  result1 = model.fitTo(*(RooDataHist*)arrayOfRooDataHist->At(0),Save(kTRUE));
result = model.fitTo(*(RooDataHist*)arrayOfRooDataHist->At(0),Extended(kTRUE),Minos(kTRUE), Hesse(kTRUE), Save(kTRUE));
}

//METHOD 2 : Standard  binned Chi2
// ---------------------------------------------------------------------
if (_scenario == 1)
{

result = model.chi2FitTo(*(RooDataHist*)arrayOfRooDataHist->At(0),Extended(kTRUE),Save(kTRUE),Minos(kTRUE),Hesse(kTRUE));
//result->Print("v");
}


//METHOD 3 : Unbinned ML with help of RooMinuit
// ---------------------------------------------------------------------
if (_scenario == 2)
{


// Construct binned likelihood
RooNLLVar nll("nll","nll",model,*(RooDataHist*)arrayOfRooDataHist->At(0),Extended(kTRUE)) ;

// Start Minuit session on NLL
RooMinuit m(nll) ;
m.migrad() ;
m.hesse() ;
result = m.save() ;
//result->Print("v") ;
}


//METHOD 4 : Binned Chi2 with help of RooMinuit
// ---------------------------------------------------------------------

if (_scenario == 3)
{


// Construct binned likelihood

RooChi2Var chi2("chi2","chi2",model,*(RooDataHist*)arrayOfRooDataHist->At(0),Extended(kTRUE),DataError(RooAbsData::SumW2)) ;


// Start Minuit session on Chi2
RooMinuit m2(chi2) ;
m2.setStrategy(2);
m2.setErrorLevel(1.);

//m2.setPrintLevel(0);
//m2.setVerbose(kTRUE);

m2.seek();
#ifdef doHESSE
m2.hesse();
#endif

result = m2.save() ;

double __tot=0;

#ifdef doMIGRAD
m2.migrad();

cout<<"\n\n +MYOUTPUT result from MIGRAD"<<"\n\n";
if ( result) result->Print("v") ;

for (int i=1; i<arrayOfNames->GetSize();i++) __tot+=(RooRealVar*)arrayOfYields->At(i))->getVal();

cout<<"\n\n VALS : "<<"\n\n";
for (int i=1; i<arrayOfNames->GetSize();i++)
{
        TString _name = ( (TObjString *) arrayOfNames->At(i))->GetString();
        double __yield= (RooRealVar*)arrayOfYields->At(i))->getVal();
        double __err = (RooRealVar*)arrayOfYields->At(i))->getError();
        cout<<_name <<" = "<<__yield<<" +- "<<__err<<"\n";
}


cout<<"\n\n RATIOS : "<<"\n\n";
for (int i=1; i<arrayOfNames->GetSize();i++)
{
        TString _name = ( (TObjString *) arrayOfNames->At(i))->GetString();
        double __yield= (RooRealVar*)arrayOfYields->At(i))->getVal();
        double __err = (RooRealVar*)arrayOfYields->At(i))->getError();
        cout<<_name <<" = "<<__yield/__tot<<" +- "<<__err/__tot<<"\n";
}


cout<<"\n\n -MYOUTPUT\n\n";

#endif

#ifdef doHESSE
m2.hesse();
#endif


#ifdef doMINOS
m2.minos();

cout<<"\n\n +MYOUTPUT result from MINOS"<<"\n\n";

result = m2.save() ;

result->Print("v") ;

__tot=0;
for (int i=1; i<arrayOfNames->GetSize();i++) __tot+=(RooRealVar*)arrayOfYields->At(i))->getVal();

cout<<"\n\n VALS : "<<"\n\n";
for (int i=1; i<arrayOfNames->GetSize();i++)
{
        TString _name = ( (TObjString *) arrayOfNames->At(i))->GetString();
        double __yield= (RooRealVar*)arrayOfYields->At(i))->getVal();
        double __err = (RooRealVar*)arrayOfYields->At(i))->getError();
        cout<<_name <<" = "<<__yield<<" +- "<<__err<<"\n";
}


cout<<"\n\n RATIOS : "<<"\n\n";
for (int i=1; i<arrayOfNames->GetSize();i++)
{
        TString _name = ( (TObjString *) arrayOfNames->At(i))->GetString();
        double __yield= (RooRealVar*)arrayOfYields->At(i))->getVal();
        double __err = (RooRealVar*)arrayOfYields->At(i))->getError();
        cout<<_name <<" = "<<__yield/__tot<<" +- "<<__err/__tot<<"\n";
}




cout<<"\n\n -MYOUTPUT \n\n";
#endif

cout<<"Optimal Chi2 from test "<<chi2.getVal()<<endl;
//result = m2.fit("svh");
//result->Print("v") ;

if (result)
{
 double dataFitNll = result->minNll() ;
cout<<"\n\n +MYOUTPUT \n\n";
	cout<<"Chi2 minimum "<<dataFitNll<<"\n";


}

///Let's do scan!
/*
Double_t _vals[40];
Double_t _errs[40];


for (int ii=1; ii<arrayOfNames->GetSize();ii++)
{
	_vals[ii-1]=(RooRealVar*)arrayOfYields->At(ii)->getVal();
	_errs[ii-1]=(RooRealVar*)arrayOfYields->At(ii)->getError();
	


}
*/

for (int ii=1; ii<arrayOfNames->GetSize();ii++)
{
int numVar=ii;					

	TString _name = ( (TObjString *) arrayOfNames->At(numVar))->GetString();                
	cout<<"\n"<<"Chi2 scan "<<_name<<"\n";


Double_t el1,el2,yel1,yel2,xmin,ymin;


TGraph * gr8 = RooFitChi2Analysis(&chi2, (RooRealVar*)arrayOfYields->At(numVar),3,0,dataYield,el1,el2,yel1,yel2,xmin,ymin);
if (gr8) {
gr8->SetName(Form("Chi2_%s",_name.Data()) );
gr8->SetTitle(Form("Chi2_%s",_name.Data()) );
arrayOfCanvas->Add (new TCanvas (Form("Chi2Scan_%s",_name.Data()),Form("Chi2Scan_%s",_name.Data())) );
TGraphLinePlot(gr8,kBlue, el1,el2,xmin);
}

cout<<"at optimal parameters: "<<"\n";
cout<<"Estimated el1 = "<<el1<<"  el2 = "<<el2<<"\n";
cout<<"Estimated yel1 = "<<yel1<<"  yel2 = "<<yel2<<"\n";
cout<<"Estimated xmin = "<<xmin<<"  ymin = "<<ymin<<"\n";
cout<<"Estimated xmin -el1 = "<<xmin-el1<<"\n";
cout<<"Estimated el2 - xmin = "<<el2-xmin<<"\n";

TGraph * gr9 = RooFitMinuitChi2Analysis( &m2,&chi2, (RooRealVar*)arrayOfYields->At(numVar),arrayOfYields,3,0,dataYield,el1,el2,yel1,yel2,xmin,ymin);

if (gr9) {
gr9->SetName(Form("Chi2_2_%s",_name.Data()) );
gr9->SetTitle(Form("Chi2_2_%s",_name.Data()) );

gr9->SetMarkerColor(kRed);
gr9->SetMarkerSize(1.1);
gr9->SetMarkerStyle(22);

//arrayOfCanvas->Add (new TCanvas (Form("Chi2Scan2_%s",_name.Data()),Form("Chi2Scan_%s",_name.Data())) );
TCanvas * ccc = (TCanvas*)arrayOfCanvas->FindObject(Form("Chi2Scan_%s",_name.Data()));
if (ccc) ccc->cd();
if (gr8) TGraphLinePlot(gr9,kRed, el1,el2,xmin,gr8->GetYaxis(),"P same",0.4);
}

cout<<"ala Minos: "<<"\n";
cout<<"Estimated el1 = "<<el1<<"  el2 = "<<el2<<"\n";
cout<<"Estimated yel1 = "<<yel1<<"  yel2 = "<<yel2<<"\n";
cout<<"Estimated xmin = "<<xmin<<"  ymin = "<<ymin<<"\n";
cout<<"Estimated xmin -el1 = "<<xmin-el1<<"\n";
cout<<"Estimated el2 - xmin = "<<el2-xmin<<"\n";



} ///for

cout<<"\n\n -MYOUTPUT \n\n";


///Lets's do contour plots

for (int ii=1; ii<arrayOfNames->GetSize();ii++)
for (int jj=ii+1; jj<arrayOfNames->GetSize();jj++)
{

        TString _name1 = ( (TObjString *) arrayOfNames->At(ii))->GetString();
        TString _name2 = ( (TObjString *) arrayOfNames->At(jj))->GetString();
        cout<<"\n"<<"Contour "<<_name1<<"_"<<_name2<<"\n";

RooRealVar * var1 = (RooRealVar*)arrayOfYields->At(ii);
RooRealVar * var2 = (RooRealVar*)arrayOfYields->At(jj);
// Make contour plot of mx vs sx at 1 sigma
  RooPlot* frame = m2.contour(*var1,*var2,1,-1) ;

TGraph * gr10=0;
for (int j=0; j<frame->numItems(); j++) {
    TString obj_name=frame->nameOf(j);
        if (obj_name=="") continue;
    TObject *obj = frame->findObject(obj_name);
    if (obj!=0 && obj->IsA()->InheritsFrom(TGraph::Class()) ) {
	//	cout<<"ob-> "<<obj_name<<"\n";
	//	cout<<"Class->"<<obj->IsA()->GetName();
		gr10= obj;

        }
}

if (gr10) {
arrayOfCanvas->Add (new TCanvas (Form("CONT0_%s_%s",_name1.Data(),_name2.Data()), Form("CONT0_%s_%s",_name1.Data(),_name2.Data()) ));
frame->Draw();
}

arrayOfCanvas->Add (new TCanvas (Form("CONT_%s_%s",_name1.Data(),_name2.Data()), Form("CONT_%s_%s",_name1.Data(),_name2.Data()) ));
gr10 = RooFitCreateContour(&m2,&chi2,var1,var2,arrayOfYields,2,1,0,dataYield);
if (gr10) {
arrayOfCanvas->Add (new TCanvas (Form("CONT2_%s_%s",_name1.Data(),_name2.Data()), Form("CONT2_%s_%s",_name1.Data(),_name2.Data()) ));
TGraphEllipsePlot(gr10,kBlue,kRed,var1->getVal(),var2->getVal());
//gr10->Draw("a*");

} //if




} ///for



}

if (result!=0) 
{

result->Print("v") ;

///Correlations  

  gStyle->SetOptStat(0) ;
  gStyle->SetPalette(1) ;
  TH2* hcorr = result->correlationHist() ;

///For plotting Correlation
arrayOfCanvas->Add (new TCanvas(TString("correlation")+TString("_canva"),TString("correlation")+TString("_canva")));
((TCanvas*)arrayOfCanvas->FindObject(TString("correlation")+TString("_canva")))->cd();

((TCanvas*)arrayOfCanvas->FindObject(TString("correlation")+TString("_canva")))->SetLeftMargin(10);
hcorr->GetYaxis()->SetTitleOffset(1.4) ;
hcorr->GetYaxis()->SetLabelSize(0.025) ;
hcorr->GetYaxis()->SetLabelOffset(-0.025) ;
//hcorr->GetXaxis()->SetLabelSize(0.025) ;
hcorr->Draw("colz");

}

///Calculate final Chi2 and prob

cout<<endl<<endl<<endl;

RooChi2Var chi2_check("chi2_check","chi2",model,*(RooDataHist*)arrayOfRooDataHist->At(0),Extended(kTRUE),DataError(RooAbsData::SumW2)) ;
cout<<std::scientific <<" 1.1) Chi2 value = " <<chi2_check.getVal() << endl ;
cout<<std::scientific <<" 1.2) ndof = " << bins-yields->getSize() << endl ;
cout<<std::scientific <<" 1.3) Chi2/ndof = " << chi2_check.getVal()/(bins-yields->getSize()) << endl ;
cout<<std::scientific <<" 1.4) Probability of Chi2 = "<<1e0- TMath::Prob(chi2_check.getVal(),bins-yields->getSize())<<endl;



//Plot results : data + templates
// ---------------------------------------------------------------------


///Create Legend
//TLegend * leg = new TLegend(0.4,0.6,0.89,0.89);
TLegend *  leg = (TLegend*) arrayOfLegends->At(arrayOfLegends->GetSize()-1); ///last;
 



RooPlot * frame = Mjj.frame() ;
((RooDataHist*)arrayOfRooDataHist->At(0))->plotOn(frame) ;
 //model.plotOn(frame,Normalization(1.0,RooAbsReal::RelativeExpected),DataError(RooAbsData::SumW2)) ;
  model.plotOn(frame,Normalization(1.0,RooAbsReal::RelativeExpected)) ;
cout<<"Number of bins in the frame "<<frame->GetNbinsX()<<endl;

///Check chi2  again!

cout<<endl<<endl<<endl;

///I don't know what frame->chiSquare() returns! People from RooTalk say that it's Chi2/ndof, but
/// RooFit developers state it as Chi2 (look at rf109_chi2residpull.C)

//double chi2_check2 = frame->chiSquare(yields->getSize());
double chi2_check2 = frame->chiSquare();
cout<<std::scientific <<"2.1) chiSquare()  = " <<chi2_check2 << endl ;
cout<<std::scientific <<"2.2) chiSquare()*ndof  = " <<chi2_check2*(frame->GetNbinsX()-yields->getSize()) << endl ;
cout<<std::scientific <<"2.3) Probability of Chi2( chiSquare)  = "<<1e0-TMath::Prob(chi2_check2,frame->GetNbinsX()-yields->getSize())<<endl;
cout<<std::scientific <<"2.4) Probability of Chi2( chiSquare*ndof)  = "<<1e0-TMath::Prob(chi2_check2*(frame->GetNbinsX()-yields->getSize()-1),frame->GetNbinsX()-yields->getSize())<<endl;


////My Chi2 calculation:

double _vals[40];
double _errs[40];


if (arrayOfHistos && arrayOfHistos->GetSize()>1  )
{


_mytemplates = new TList();
_myyields = new TArrayD(arrayOfNames->GetSize());
for (int i=0; i<arrayOfNames->GetSize();i++)
{
	TH1 * tempHis = (TH1*)arrayOfHistos->At(i);
	TString _name = ( (TObjString *) arrayOfNames->At(i))->GetString();

	double _yield=-1e10;
	
if (tempHis!=0) 
{ 
	if ( i==0) {
		_mytemplates->Add(tempHis);  ///data
		_yield=0;
		}
	else 
		{
	TH1F* _histPdf=(TH1F*) ((RooHistPdf*)arrayOfRooHistPdf->At(i))->createHistogram(_name+TString("_histpdf_validation2"),Mjj));
		_mytemplates->Add(_histPdf);	
		_yield= (RooRealVar*)arrayOfYields->At(i))->getVal();
		_vals[i] = _yield;
		_errs[i] = (RooRealVar*)arrayOfYields->At(i))->getError();
		}
	

		cout<<"I'm adding "<<((TObjString*) arrayOfNames->At(i))->GetString()<<" with yield "<<_yield<<endl;
		_myyields->AddAt(_yield,i);

} /// if tempHis

} /// for

///Plotting the ratios


double _tot=0;
for (kk=0;kk<_myyields->GetSize();kk++) _tot+=_myyields->At(kk);
 cout<<"DATA = "<<dataYield<<endl;
 cout<<"MC = "<<_tot<<endl;

for (int i=1; i<arrayOfNames->GetSize();i++)
{
        TString _name = ( (TObjString *) arrayOfNames->At(i))->GetString();
	cout<<_name<<"  ---> " <<" yield : "<< _vals[i]<<" +- "<< _errs[i]<<" ratio : "<<_vals[i]/_tot<<" +- "<<_errs[i]/_tot<<"\n";

}


 
 





//double _mychi2 = chi2test(_mytemplates,_myyields,"UU P CHI2");
double _mychi2 = chi2test2(_mytemplates,_myyields);


cout<<std::scientific <<"3.1) mychi2  = " <<_mychi2 << endl ;
cout<<std::scientific <<"3.2) mychi2/ndof  = " <<_mychi2/(((TH1*)arrayOfHistos->At(0))->GetNbinsX()-_myyields->GetSize()) << endl ;
cout<<std::scientific <<"3.3) Probability of mychi2  = "<<1e0-TMath::Prob(_mychi2, ((TH1*)arrayOfHistos->At(0))->GetNbinsX()-_myyields->GetSize())<<endl;

///Clean up
for (int i=1; i<_mytemplates->GetSize();i++)  if (_mytemplates->At(i)) delete _mytemplates->At(i);
delete _mytemplates;
delete _myyields;


} ///if arrayOfHistos

//for (unsigned int i=0;i<sizeof(histonames)/sizeof(char *);i++)
for (int i=0; i<arrayOfNames->GetSize();i++)
{

if (i==0) continue; /// skip data

TString _name = ( (TObjString *) arrayOfNames->At(i))->GetString();

if (i==0)  continue; /// skip data
 model.plotOn(frame, Components(_name+TString("_histpdf")),LineColor(colors[i]),Normalization(1.0,RooAbsReal::RelativeExpected));
}




for (int j=0; j<frame->numItems(); j++) {
    TString obj_name=frame->nameOf(j); 
	if (obj_name=="") continue;
    TObject *obj = frame->findObject(obj_name);
    if (obj!=0) {
		
	   leg->AddEntry(obj,obj_name,"l");
	}
}



///Make Sorting

TGraph _comb;
TList * lstTH1 = new TList();
for (int i=1; i<arrayOfNames->GetSize();i++)
{
TString _name = ( (TObjString *) arrayOfNames->At(i))->GetString();


TH1F* tempHis2=(TH1F*) ((RooHistPdf*)arrayOfRooHistPdf->At(i))->createHistogram(_name+TString("_"),Mjj));

if (!tempHis2) continue;
tempHis2->Scale(_vals[i]);

lstTH1->Add(tempHis2);
_comb.SetPoint(i,i,tempHis2->GetMaximum());

}



_comb.Sort(&TGraph::CompareY,kFALSE); ///descending on Y
Double_t * templPos = _comb.GetX();

arrayOfCanvas->Add (new TCanvas ("final_stack","final stack"));

TCanvas * c5 = (TCanvas*) arrayOfCanvas->Last();
TLegend * leg2 = new TLegend(0.4,0.6,0.89,0.89);
THStack *hs = new THStack("hs","final stack");
if (dataHist) leg2->AddEntry(dataHist,"data","l");

for (int i=0; i<_comb.GetN();i++)
{
if (templPos[i] == 0 ) continue; ///artifact --> TGraph always adds (0,0) point if it doesn't exist
TString _name = ( (TObjString *) arrayOfNames->At(templPos[i]))->GetString();

TH1F * tempHis3 = (TH1F*)lstTH1->At(templPos[i]-1);


if (!tempHis3) continue;


tempHis3->SetLineColor(colors[int(templPos[i])]);
tempHis3->SetFillColor(colors[int(templPos[i])]);
tempHis3->SetFillStyle(1001);

leg2->AddEntry(tempHis3,_name,"l");
hs->Add(tempHis3);

}

c5->cd();
hs->Draw("hist");
if (dataHist) {
dataHist->Draw("e0 same");
}
leg2->Draw();





((TCanvas*)arrayOfCanvas->FindObject(TString("final")+TString("_canva")))->cd();

frame->Draw();
leg->Draw();

cout<<endl<<endl<<endl;


///Print info
model.getParameters(Mjj)->Print("s") ;



cout<<endl<<endl<<endl;






// Print parameter list in LaTeX for (one column with names, one column with values)
cout<<"Latex table "<<endl<<endl<<endl;
  params->printLatex() ;


///Save figures
for (int j=0; j<arrayOfCanvas->GetEntries(); j++) {

if (j==0) continue;

TString _posfix = Form("_Btag%d_Corr%d_Scenario%d",_novl, _corr, _scenario );
((TCanvas*)arrayOfCanvas->At(j))->SaveAs(TString(arrayOfCanvas->At(j)->GetName())+_posfix+TString(".eps"));

}


///Clean

delete arrayOfHistos;
delete arrayOfRooDataHist;
delete arrayOfRooHistPdf;
//delete arrayOfCanvas;
delete arrayOfYields;
delete arrayOfLegends;
delete _combinedList;

delete yields;
delete shapes;



return;
}
